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ABSTRACT 

Very high energy (VHE, E >100 GeV) 7 -ray flaring activity of the high- 
frequency peaked BL Lac object PG 1553-1-113 has been detected by the H.E.S.S. 
telescopes. The flux of the source increased by a factor of 3 during the nights 
of 2012 April 26 and 27 with respect to the archival measurements with hint 
of intra-night variability. No counterpart of this event has been detected in 
the FermLLAT data. This pattern is consistent with VHE 7 ray flaring being 
caused by the injection of ultrarelativistic particles, emitting 7 rays at the highest 
energies. The dataset offers a unique opportunity to constrain the redshift of 
this source at z = 0.49 ± 0.04 using a novel method based on Bayesian statistics. 
The indication of intra-night variability is used to introduce a novel method to 
probe for a possible Lorentz Invariance Violation (LIV), and to set limits on the 
energy scale at which Quantum Gravity (QG) effects causing LIV may arise. 
For the subluminal case, the derived limits are Eqg i > 4.10 x 10^^ GeV and 
Eqg,2 > 2.10 X 10^° GeV for linear and quadratic LIV effects, respectively. 

Subject headings: Galaxies: active - BL Lacertae objects: Individual: PG 1553-1-113 
- Gamma rays: observations - Quantum Gravity - Lorentz invariance breaking 


1. Introduction 


Blazars are activ e galactic nuclei (AGN) with their jets closely aligned with the line 
of sight to the Earth fjUrrv fc Padovanilll995l) . Among their particularities is flux variabil- 
i ty at all wavelengths on various time scal es, from years down to (in some cases) minutes 
(IGaidos et al.lll996l : lAharonian et al.ll2007a|l . Flaring activity of blazars is of great interest for 
probing the source-intrinsic physics of relativistic jets, relativistic particle acceleration and 
generation of high-energy radiation, as well as for conducting fundamental physics tests. On 
the one hand, exploring possible spectral variability between flaring and stationary states 
helps to understand the electromagnetic emission mechanisms at play in the jet. On the 
other hand, measuring the possible correlation between photon energies and arrival times 
allows one to test for possible Lorentz invariance violation (LIV) leading to photon-energy- 
dependent variations in the speed of light in vacuum. 


Holesovickach 2, 180 00 Prague 8, Czech Republic 

^"'’GRAPPA, Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH 
Amsterdam, The Netherlands 
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Located in the Serpens Caput constellation, PG 1553+113 was discovered bv iGreen et al 


fll986l) . who hrst classihed it as a BL Lac ob ject. Later the class ihcation was refined to a 


high-frequency peaked BL Lac object (HBL, iGiommi et ahl 1199511 . PG 1553 +113 exhibits 


a high X-ray to radio flux (log(F 2 ohz) > —4.5. lOsterman et ahl 1200611 . which places 


it among the most extreme HBLs (IRector et ahl 120031) . The object was observed in X 


0.3 X 10 


-11 


erg cm ^ s ^ 


Osterman et al. 

2006 

) to 3.5 X 10 ^^ergcm ^s ^ ( 

Reimer et al. 


20081 1 but no fast variability (in the sub-hour time scale) has been detected so far. 


PG 1553+113 was discove red at very high energies (VHE, E >100 GeV) by H.E.S.S 
flAharonian et ahl l2006ai 1200811 with a photon index of P = 4.0 ± 0.6. At h igh energies 


(HE, 100 MeV < E <300 GeV) the source has been detected by the Ferm 2 -LAT flAbdo et al. 


2009bl. l2010a|l with a very hard photon index of P = 1.68 ± 0.03, making this object the 
one with the largest HE ~ VHE spectra l break (AP ^ 2.3 ) ever measured. No variability 
in Fermz-LAT was found by lAbdo et ahl fl2009bl l2010alj on daily or weekly time scales, but 
using an extended data set of 17 months. lAleksic et al. ( 2012 ) reported variability above 
1 GeV with flux variations of a factor of ~ 5 on a yearly time scale. 


With 5 years of monitoring data of the MAGIG telescopes, lAleksic et ahl fj 2012 l) discov¬ 
ered variability in VHE 7 rays with only modest flux variations (from 4 to 11 % of the Grab 
Nebula flux). In addition to the high X-ray variabil ity, this behavior can be interpreted as 
evidence for Klein-Nishina effects flAbdo et al.ll2010al) in the framework of a s ynchrotron self- 
Gompton model. The sou rce underwent VHE 7 -ray flares in 2012 March flGortinall2012al) 
and April (IGortinal l2012bll . detected by the MAGIC telescopes. During the March flare, 
the source was at a flux level of about 15% of that of the Crab Nebula, while in April it 
reached ~ 50%. During those VHE 7 -ray flares, also a brightening in X-ray, UV and op¬ 
tical wavelengths has been noticed by the MAGIC collabo ration. A detailed study of the 
MAGIC telescopes and multi-wavelength data is in press (lAleksic et ahl 1201411 . The latter 
event triggered the H.E.S.S. observations reported in this work. 

Despite several attempts to measure it, the redshift of PG 1553+113 still suffers from un- _ 

certa inties. Different attempts, including optical spectroscopy ( Treves et ahl 2007 : Aharonian et al. 


20081) or comparison s of the HE and VHE spectra of PG 1553+113 (iPrandini et al 


Sanchez et al. 


2009 


201311 . were made. Based on t he assumption t hat t he EBL-corrected VHE 


spectral index is equal to the Fermz-LAT one, iPrandini et ahl (1200911 derived an upper limit 
oi z < 0.67. Gomparing PG 1553+113 statistically with other known VHE emitters and tak- 
ing into account a p ossible intrinsic 7 -ray spectral break through a simple emission model, 
Sanchez et ahl (l2013ll constraine d the r edshift to be below 0.64. The best estimate to-date 
was obtained bv iDanforth et ahl (1201011 who found the redshift to be between 0.43 and 0.58 
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using far-ultraviolet spectroscopy. 

This paper concentrates on the HE and VHE emission of PG 1553-1-113 and is divided 
as follow: Sections 12.11 and 12.21 present the H.E.S.S. and Ferm 2 -LAT analyses. The discus¬ 
sion, in section |3l includes the determination of the redshift using a novel method and the 
constraints derived on LIV using a modified likelihood formulation. Throughout this paper 
a ACDM cosm ology with Hn = 70.4 ± 1.4 kms“^ Mpc“^, Qm = 0.27±0.03, Ha = 0.73 ±0.03 
from WMAP (iKomatsu et ahl 1201 ill is assumed. 


2. Data analysis 


2.1. H.E.S.S. observations and analysis 

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes located in the 
Khon ias highland in Namibia ('23°16H8^^ S, 16°3T 01" E), at an altitude of 1800 m above sea 
level ( Hinton fc the H.E.S.S. Collaboration 2004). The fifth H.E.S.S. telescope was added 
to the system in 2012 July and is not used in this work, reporting only on observations prior 
to that time. 


PG 1553-1-113 was observed with H.E.S.S. in 2005 and 2006 ( Aharonian et al.ll2008 ). No 
variability was found in these observations, which will be referred to as the “pre-flare” data 
set in the following. New observations were carried out in 2012 Apr il after flaring activity 
at VHE was reported by the MAGIC collaboration (“flare” data set. ICortinall2012bl) . 


_ The pre-flare data set is composed of 26.4 live time hours of good-quality data flAharonian et al. 

2006b|). For the flare period, eight runs of ~ 28 minutes each were taken during the nights 


of 2012 April 26 and 27, corresponding to 3.5 hours of live time. All the data were taken in 
wobble mode, for which the source is observed with an offset of 0?5 with respect to the center 
of the instrument’s held of view yielding an acceptance-corrected live time of 24.7 hours and 
3.2 hours for the pre-hare and hare data sets, respectively. 


Data were analyzed using the Model analysis fide Naurois fc Rollandll2009l) with Loose 
cuts. This method-based on the comparison of detected shower images with a pre-calculated 
model-achieves a better rejection of hadronic air showers and a better sensitivity at lower 
energies than analysis methods based on Hillas parameters. The chosen cuts, best suited 
for sources with steep spectra such as PG 1553-|-lHii|, require a minimum image charge of 
40 photoelectrons, which provides an energy threshold of ~ 217 GeV for the pre-hare and 


^PG 1553-1-113 has one of the steepest spectra measured at VHE. 
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240 GeV for the flare data set@. All the results presented in thi s pape r were cross-checked 
with the independent analysis chain described in iBecherini et al.l (120 111) . 


Events in a circular region (ON region) centered on the r adio position of the source, 
'^j 2 ooo = 15*^55™43.04®, ( 5 j 2 ooo = 11°11'24.4" f Green et al. 1986 ). with a maximum squared 
angular distance of 0.0125 deg^, are used for the analysis. In order t o estimate the background 
in this region, the reflected background method (IBerge et al.ll2007l ) is used to dehne the OFF 
regio ns. The excess of 7 rays in the ON region is statistically highly signihcant (ILi fc Ma 
19831 ): 21 . 5(7 for the pre-flare period and 22.0a for the flare. Statistics are summarized in 
Table [H 


The differential energ y spectrum of the VHE 7 -ray emission has been derived using 
a forward-folding method fjPiron et al.ll200ll) . For the observations prior to 2012 April, a 
power law (PWL) model htted to the data gives a of 51.7 for 40 degrees of freedom (d.o.f., 
corresponding to a probability of P ^2 = 0.10). The values of the spectral parameters 
(see Table El) are comp atible with previous analyses by H.E.S.S. covering the same period 
( Aharonian et al. 20081) . A log-parabola (LP) models, with a x^ of 37.5 for 39 d.o.f. (P ^2 = 
0.54), is found to be preferred over the PWL model at a level of 4.3 a using the log-likelihood 
ratio test. Note that sys tematic uncertainties, presented in Table [21 have been evaluated by 


Aharonian et al.l fl2006bl) for the PWL model and using the jack-knife method for the LP 
model. The jack-knife method consist in removing one run and redoing the analysis. This 
process is repeated for all runs. 


For the flare data set, the log-parabola model does not signihcantly improve the ht and 
the simple PWL model describes the data well, with a x^ of 33-0 for 23 d.o.f. {P ^2 = 0.08). 
Table [2] contains the integral fluxes above the reference energy of 300 GeV. The flux increased 
by a factor of ~ 3 in the flare data set compared to the pre-flare one with no sign of spectral 
variations (when comparing power law hts for both data sets). The derived spectra and error 
contours for each data set are presented in Fig. [H where the spectral points obtained from 
the cross-check analysis are also plotted. 

To compute the light curves, the integrated flux above 300 GeV for each observation 
run was extracted using the corresponding (pre-flare or flare) best ht spectral model. A ht 
with a constant of the run-wise light curve of the entire (pre-hare+hare) data set, weighted 
by the statistical errors yields a x^ of 123.2 with 68 d.o.f. (P ^2 = 6.6 x 10“®). Restricting the 


^The difference of energy threshold between the two data set is due to the changing observation conditions, 
e.g., zenith angle and optical efficiency. 


^The log-parabola is defined by dN/dE — (E/Eq) “ biog(E/Eo) 
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Table 1. Summary of the statistics for both data sets (hrst column). The second and 
third columns give the number of ON and OFF events. The 4th column gives the ratio 
between ON and OFF exposures (r). The excess and the corresponding signihcance are 
given, as well as the energy threshold and the mean zenith angle of the source during the 
observations. The last column presents the probability of the flux to be constant within the 

observations (see text). 


Data set 

ON 

OFF 

r 

Excess 

Significance 

^th [GeV] 

Zenith angle 

pest 

Pre-Flare 

2205 

13033 

0.100 

901.7 

21.5 

217 

34° 

0.77 

Flare 

559 

1593 

0.105 

391.2 

22.0 

240 

52° 

3.3 X 10-3 


True energy (TeV) 


True energy (TeV) 



Reconstructed energy (TeV) 


Reconstructed energy (TeV) 


Fig. 1.— Differential fluxes of PG 1553+113 during the pre-flare (left) and flare (right) 
periods. Error contours indicate the 68 % uncertainty on the spectrum. Uncertainties on 
the spectral points (in black) are given at 1 cr level, and upper limits are computed at the 
99 % conhdence level. The gray squares were obtained by the cross-check analysis chain and 
are presented to visualize the match between both analyses. The gray error contour on the 
left panel is the best-£t power law model. The lower panels show the residuals of the £t, 
i.e. the difference between the measured (nobs) and expected numbers of photons (nmodei), 
divided by the statistical error on the measured number of photons 
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Fig. 2.— H.E.S.S. light curve of PG 1553+113 during the 2 nights of the flare period. The 
continuous line is the measured flux during the flare period while the dashed one corresponds 
to the pre-flare period (see Table [2] for the flux values). Gray areas are the 1 a errors. 


analysis to the pre-flare data set only, the £t yields a of 51.76 with 60 d.o.f. (P ^2 = 0.77), 
indicating a gain a f l ux incr ease detected by H.E.S.S. at the time of the flaring activity 
reported by iGortinal f)2012bl) . 


Figured shows the light curve during the flare together with the averaged integral fluxes 
above 300 GeV of both data sets. A £t with a constant to the H.E.S.S. light curve during 
the hrst night yields a x^ of 20.76 for 6 d.o.f. (P ^2 = 2.0 x 10 ^), indicati ng intra-night 


variability. This is also supported by the use of a Bayesian block algorithm flScargldllOOSll 
that hnds three blocks for the 2 nights at a 95% conhdence level. 
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2.2. Fermz-LAT analysis 


The Fermi Lar ge Area Telescope (LAT) is detector converting 7 ray to e’''e“ pairs 
( Atwood et al. 2009 '). The LAT is sensitive to 7 rays from 20 MeV to > 300 GeV. In snrvey 
mode, in which the bnlk of the observations are performed, each sonrce is seen every 3 honrs 
for approximately 30 minntes. 


The Ferm«-LAT data and software are available from the Fermi Science Snpport Cen- 
tejjl. In this work, t he Science Tools V9R32P5 were nse d with the Pass 7 reproc essed data 


(IBregeon et ahll2013f) . specihcally SOURCE class event flAckermann et alJ l2012al) . with the 
associated P7REP_S0URCE_V15 instrnment response fnnctions (IRFs). Events with energies 
from 300 MeV to 300 GeV were selected. Additional cnts on the zenith angle (< 100°) and 
rocking angle (< 52°) were applied as recommended by the LAT collaboratioE@ to reduce 
the contamination from the Earth atmospheric secondary radiation. 


The analysis of the LAT data was performed using the Enrico Python package fjSanchez fc Dei! 


2 OI 3 I ). The sky model was dehned as a region of interest (ROI) of 15° radius with PG 1553+113 
in the center and additional point-like sources from the internal 4-years source list. Only the 
sources within a 3° radius around PG 1553+113 and bright sources (integral flux greater that 


5 X 10“'phcm“^ s“^) had their parameters free to vary during the likelihood minimization. 
The template hies isotrop_4years_P7_V15_repro_v2_source.txt for the isotropic diffuse 
component, and template_4vears_P7_vl5_ repro_v2. f its for t he standard Galactic model, 
were included. A binned likelihood analysis flMattox et al.lIlQQGI) . implemented in the gtlike 
tool, was used to hnd the best-ht parameters. 


As for the H.E.S.S. data analysis, two spectral models were used: a simple PWL and a 
LP. A likelihood ratio test was used to decide which model best describes the data. Table [3] 
gives the results for the two time periods considered in this work, and Figure [3] presents 
the 7 -ray SEDs. The hrst one (pre-hare), before the H.E.S.S. exposures in 2012, includes 
more than 3.5 years of data (from 2008 August 4 to 2012 March 1). The best ht model is 
found to be the LP (with a T^^ of 11.3, ~ 3.4cr). The second period (hare) is centered on 
the H.E.S.S. observations windows and lasts for seven days. The best ht model is a power 
law, the hux being consistent with the one measured during the hrst 3.5 years. Data points 
or light curves were computed within a restricted energy range or time range using a PWL 


"^http: //f ermi . gsfc.nasa.gov/ssc/data/analysis/ 

“http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/index.html 

®Here the TS is 2 times the difference between the log-likelihood of the fit with a LP minus the log- 
likelihood with a PL. 
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model with the spectral index frozen to 1.70. 


To precisely probe the variability in HE 7 rays, seven-day time bins were used to compute 
the light curve of PG 1553-1-113 in an extended time window (from 2008 August 4 to 2012 
October 30), to probe any possible delay of a HE flare with respect to the VHE one. While 
the flux of PG 1553-1-113 above 300 M eV is found to be vari able in the whole period with a 
variability index of F^ar = 0.16 ± 0.04 fjVaughan et al.ll2003l) . there is no sign of any flaring 
activity around the 2012 H.E.S.S. observations. This result has been conhrmed by using 
the Bayesian block algorithm, which Ends no block around the H.E.S.S. exposures in 2012. 
Similar results were obtained when considering only photons with an energy greater than 
1 GeV. No sign of enhancement of the HE flux associated to the VHE event reported here 
was found. This might be due to the lack of statistic at high energy in the LAT energy 
range. 



Fig. 3.— Spectral energy distribution of PG 1553-1-113 in 7 rays as measured by the Fermi- 
LAT and H.E.S.S. Red (blue) points and butterflies have been obtained during the flare 
(pre-flare) period. The Fermi and H.E.S.S. data for the pre-flare are not contemporaneous. 
H.E.S.S. data were taken in 2005-2006 while the Fermi data were taken between 2008 and 
2012 . 
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Table 2. Summary of the fitted spectral parameters for the pre-flare and the flare data 
sets and the corresponding integral flux / calculated above 300 GeV. The last column gives 

the decorrelation energy. 


Data Set (Model) 


Spectral Parameters 

/ (E>300 GeV) 
[10“^^ ph cm“^ s“^] 

Edec 

[GeV] 

Pre-Flare (PWL) 

P 

= 4.8 ± 0.2stat ± 0.2sys 

4.4 ± 0.4stat ± 0.9sys 

306 

Pre-Flare (LP) 

a 

b 

= 5.4 ± 0.4stat ± O.lsys 
= 4.0 ± 1.4<,tat ± 0.2gys 

5.0 ± O.Ostat ± l.Osys 


Flare (PWL) 

F 

= 4.9 ± 0.3stat ± 0.2sys 

15.1 ± 1.3stat ± 3.0sys 

327 


Table 3. Results of the Fermi-LAT data analysis for the pre-flare and flare periods. For 
the latter, the analysis has been performed in two energy ranges (see 13.2p . The hrst 
columns give the time and energy windows and the third the corresponding test statistic 
(TS) value. The model parameters and the flux above 300 MeV are given in the last 


columns. The systematic uncertain 


les were computed using the IRFs bracketing method 


flAbdo et al.ll2009al) . 


MJD range 

Energy range 
[GeV] 

TS 


Spectral Parameters 

/(E> 300MeV) 
10“®[ph cm“^ s“^] 

54682-55987 

0.3-300 

7793.7 

a = 

= 1.49 ± 0.06stat ± O.Olgys 

2.82 ± O.lstat ± 0.2sys 




b = 

= 3.8 ± l.lstat ± O.lsys 


56040-56047 

0.3-300 

43.8 

F = 

= 1.78 ± 0.24stat ± O.Olsys 

3.5 ± 1.3stat ± 0.3sys 

56040-56047 

0.3-80 

44.5 

F = 

= 1.72 ± 0.26stat ± O.Olsys 

3.4 ± 1.3stat ± 0.3sys 
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Discussion of the results 


3.1. Variability in 7 -rays 

The VHE data do not show any sign of variation of the spectral index (when comparing 
flare and pre-flare data sets with the same spectral model), and in HE no counterpart of this 
event can be found. The indication for intra-night variability is similar to other TeV HBLs 
(Mrk 421, Mrk 501 or PKS 2155-304) with, in this case, flux variations of a factor 3. 

As notic ed in previous work s, PG 1553-1-113 presents a sharp break between the HE and 
VHE ranges (lAbdo et al.ll2010a|) and the peak position of the 7 -ray spectrum in the 
representation is located around 100 GeV. This is confirmed by the fact that the log-parabola 
model better represents the pre-flare period in HE. Nonetheless, the precise location of this 
peak cannot be determined with the FermhLAT data only. Gombining both energy ranges 
and fitting the HE and VHE data points with a power law with an exponential cutoflLl allows 
us to determine the peak position for both time periods. The functional form of the 

model is 

E 


r.2dN 

-= N 

dE 


100 GeV 


-r 


exp{-E/Ec). 


For this purpose, Fermi- LAT and H. E. S. S . syste matic uncertainties were taken into 
account in a similar way as in lAbramowski et al.l (120141 1 and added quadr atically to the sta¬ 


tistical errors. The Fermi-LAT systematic uncertainties were estimated bv I Acker maun et ah 
(2012a) to be 10 % of the effective area at 100 MeV, 5 % at 316 MeV and 15 % at 1 TeV and 


above. For the VHE 7 -ray range, they were taken into account by shifting the energy by 10 %. 
This effect translates into a systematic uncertainty for a single point of cr(f)sys = 0.1 ■ df/dE 
where / is the differential flux at energy E. 

The results of this parameterization are given in Table 01 Using the pre-flare period, 
the peak position is found to be located at log^g(Emax/l GeV) = 1.7 ± 0.2stat ± 0.4sys with 
no evidence of variation during the flare and no spectral variation. This is consistent with 
the fact that no variability in HE 7 rays was found during the H.E.S.S. observations. This 
is also in ag reement with the f act that HBLs are less variable in HE 7 rays than other BL 
Lac objects flAbdo et al.ll2010bl ). while numerous flares have been reported in the TeV band. 


fit with a LP model has been attempted, but the power law with an exponential cutoff leads to a 
better description of the data. 















Table 4. Parametrization results of the two time periods (first column) obtained by combining H.E.S.S. and 
Fermz-LAT. The second column gives the normalization at 100 GeV, while the third and the fourth present the 
spectral index and cut-off energy of the htted power law with an exponential cut-off. The last column is the peak 

energy in a zz/(z/) representation. 


Period 

N (E=100 GeV) 
10“^^ [erg cm“^ s“^] 

F 


logio(Ee/l GeV) 

logio(Emax/l GeV) 

Pre-Flare 

Flare 

9.6 ± O.Tstat ± l-7sys 
13.0±3.5,tat±5.7,y, 

1.59 ± 0.02,tat 
1.56 ± 0.08,tat 

± 0.03sys 
± O.llsys 

2.03 ± 0.02gtat ± O.Odsys 
2.16 ± O.Odstat ± O.OOsys 

1.7 ± 0.2stat ± 0.4sys 

1.8 ± 0.7stat ± 1.3sys 
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3.2. Constraints on the redshift 


The extragalactic background light (EBL) is a held of UV to far infrared photons pro- 
duced by the thermal emission from stars and reprocessed starlight by dust in galaxies (see 
Hauser &: Dwekll200ll. for a review) that interacts with very high energy 7 rays from sources 
at cosmological distances. As a consequence, a source at redshift 2 ; exhibits an observed 
spectrum (pobsiE) = (pint{E) x where (p-mtiE) is the intrinsic source spectrum and 

r is the optical depth due to interaction with the EBL. Since the optical depth increases 
with increasing 7 -ray energy, th e integral hux is lowered a nd the spectral index is increasecl^. 
In the following, the model of iFranceschini et ahl (1200811 was used to compute the optical 
depth r as a function of redshift and energy. In this section, the data taken by both instru¬ 
ments during the flare period are used, with the Fermi-LAT analysis restricted to the range 
300 MeV< E <80 GeV (see Table [3] for the results). In the modest redshift range of VHE 
emitters detected so far {z < 0.6), the EBL absorption is negligible below 80 GeV ~ 0.1 
at 80 GeV for x = 0.6). 


A measure of t he EB L energy density was obtained by I Ackerman n et ahl fj2012bfl and 
Abramowski et al.l (j2013bl) based on the spectra of sources with a known x. In the case 
of PG 1553-1-113, for which the redshift is unknown, the effects of the EBL on the VHE 
spectrum might be used to derive constraints on its distance. Ideally, this would be done 
by comparing the observed spectrum with the intrinsic one but the latter is unknown. The 
Fermi-LAT spectrum, derived below 80 GeV, can be considered as a proxy for the intrinsic 
spectrum in the VHE regime, or at least, as a solid upper limit (assuming no hardening of 
the spectrum). 


Following the method used bv I Abramowski et al.l fj2013al) . it has been assumed that the 
intrinsic spectrum of the source in the H.E.S.S. energy range cannot be harder than the 
extrapolation of the Fermi-LAT measurement. From this, one can conclude that the optical 
depth cannot be greater than given by: 


rmax(h^) = In 




( 1 ) 


(1 - Q:)(0obs - 1.64A(/)obs) 
where 0int is the extrapolation of the Fermi-LAT measurement towards the H.E.S.S. energy 
range. 0obs ± A0obs is the measured flux by H.E.S.S. The factor (1 — a) = 0.8 accounts 
for the systematic uncertainties of the H.E.S.S. measurement and the number 1.64 has been 


®For sake of simplicity it is assumed here that the best-fit model is a power law, an assumption which is 
true for most of the cases due to limited statistics in the VHE range 





















17 


calculated to have a confidence level of 95% flAbramowski et al.ll2013al) . The comparison is 
made at the H.E.S.S. decorrelation energy where the flux is best measured. 



Fig. 4.— Values of Tmax as a function of the photon energy. The black line is the 95% UL 
obtained with the H.E.S.S. data a nd the red line is the optical depth computed with the 
model of iFranceschini et al.l (120081) for a redshift of 0.43. The blue line is the decorrelation 
energy fo the H.E.S.S. analyse. The gray lines are the value of optical depth for different 
redshift. 


Figure 0] shows the 95 % UL on Tmax- The resulting upper limit on the redshift is 
z < 0.43. This method does not allow the statistical and systematic uncertainties of the 
FermLLAT measurement to be taken in to account and does not tak e advantage of the spectral 
features of the absorbed spectrum (.see lAbramowski et al.ll2013b|) . 


A Bayesian approach has been developed with the aim of taking all the uncertainties 
into account. It also uses the fact that EBL-absorbed spectra are not strictly power laws. 
The details of the model are presented in Appendix and only the main assumptions and 
results are recalled here. Intrinsic curvature between the HE and VHE ranges that naturally 
arises due to either curvature of the emitting distribution of particles or emission effects (e.g. 
Klein-Nishina effects) is permitted by construction of the prior (Eq. lAip : A spectral index 
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softer than the Fermz-LAT measurement is allowed with a constant probability, in contrast 
with the previous calculation. It is assumed that the observed spectrum in VHE 7 rays 
cannot be harder than the Fermi-\jPC£ measurement by using a prior that follows a Gaussian 
for indices harder than the Fermi-LKI one. The prior on the index is then: 


if r < Tpermi and 


F(r) OC ^/^(r, rpermi, fr) 
F(r) oc 1 


( 2 ) 


otherwise. TFermi is the index measured by Fermi-LAT and ar is the uncertainty on this 
measurement that takes all the systematic and statistical uncertainties into account. 



Fig. 5.— Posterior probability density as a function of redsh ift (red). The blue area rep¬ 
resents the redshift r ange estimated bv | pan forth et al.l ( 120101 ) while the green dashed line 
indicates the limit of ISanchez et al.l (120131) . 


The most probable redshift f ound with this inethod is z = 0.49±0.04, in good agreement 
with the independent measure of iDanforth et ahl (120101) . who constrained the distance to be 
between 0.43 < .^ < 0.58. Figure Ogives the posterior probability obtained with the Bayesian 
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method compared with other measurements of Lower and upper limits at a conhdence 
level of 95 % can also be derived as 0.41 < z < 0.56. Note that this method allows the 
systematic uncertainties of both instruments {Fermi-LAT and H.E.S.S.) to be taken into 
account. The spectral index obtained when htting the H.E.S.S. data with an EBL absorbed 
PWL using a redshift of 0.49 is compatible with the Fermi measurement below 80 GeV. 


3.3. Lorentz Invariance Violation 


As stated in section EH the H.E.S.S. data of the flare show a indication of intra-night 
variability, which is used here to test for a possible Lorentz Invariance Violation (LIV). Some 
Quantum Gravity (QG) models predict a change of the speed of light at energies c lose t o the 
Planck s c ale 10^® GeV). A review of such models can be found in Mattingly ( 20051 ) and 


Liberate (120131 1. An energy-dependent dispersion in vacuum is searched for in the data by 


testing a correlation between arrival times of the photons and their energies. For two pho¬ 
tons with arrival times ti and t 2 and energies Ei and E 2 , the dispersion parameter of order 


n is dehned as = 


t2-t\ 


At 

A{E") • 


Here only the linear (n = 1) and quadratic (n = 2) dis¬ 


persion parameters are calculated. Assuming no intrinsic spectral variability of the source, 
the dispersion can be related to the normalized distance of the source corrected for 


the expansion of t 
expected to occur 


le Universe and an energy Eqq at which Quantum Gravity effects are 


Jacob fc Piranll2008l) : 


_ (1 + ^) .ON 

A(E^) “ ®^E^g2Ho^” 

where Hq is the Hubble constant and s± = —1 (resp. -|-1) in the superluminal (resp. sublu¬ 
minal) case, in which the high-energy photons arrive before (resp. after) low-energy photons. 
The normalized distance is calculated from the redshift of the source and the cosmo¬ 
logical parameters flm, Ga given in the introduction: 


r {l + z^dz' 

= / - , = 

Jo \/Gm(l + z')^ + Ga 

Using the central value of z = 0.49 determined in section [321 the distance Knf 01 n 
2 is Ki = 0.541 and ^2 = 0.677. 

First, the dispersion measurement method will be described. It will then be applied to 
the H.E.S.S. flare dataset (MG simulations and original dataset), in order to measure the 
dispersion and provide 95 % 1-sided lower and upper limits on the dispersion parameter r„. 
These limits on will lead to lower limits on Eqq using equation [3l 


(4) 

= 1 and 
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3.3.1. Modified maximum likelihood method 


A maximum likelihood method, 
culate the dispersion parameter Tr,. 

bllowine: Martinez & Errando f2009h was used to cal- 

Albert et ah (2008) applied this method to a flare of 

Mkn 501, while Abramowski et ah 

2011 

applied it to a flare of PKS 2155-304. More re- 

centlv. it was used bv Vasileiou et ah 

(2011 

i) to analyse Fermi data of four gamma-ray bursts. 


The data from Cherenkov telescopes is contaminated by vr® decay from proton showers, 
misidentihed electrons, or heavy elements such as helium. In the case of PG 1553+113, and 
contrary to previous analyses, this background is not negligible: the signal-over-background 
ratio S/B is about 2, c ompared to 300 for the PKS 2155-304 flare event of July 2006 
( Aharonian et ah 2007b[l . The background was included in the formulation of the proba¬ 
bility density function (PDF) used in a likelihood maximization method. Given the times 
ti and energies Ei of the gamma-like (ON) particles received by the detector, the unbinned 
likelihood, function of the dispersion parameter is: 


"■ON 


L{rn) = JJP(Ai,+|r„). 

i=l 

The PDF P{Eifii\Tn) associated with each ON event is composed of two terms: 


(5) 


-P(GiI+1) ■ T3ig(F/j, tjIT^) + (1 Wg) ■ TBi;g(F/j, tj) 


( 6 ) 


with 

PsifiEi, ti\Tn) = I ' Aes{Ei, ti) Asig(Ai) Esifiti - Tn ' Af) 

A [Tn) 

PskgiEi, ti) = Agff (Aj, ti) A^^fiEi) Agijg 

noN — « ^OFF 
Ws = -. 

The PDF Pgig includes the emission time distribution of the photons -^Sig determined from 
a parametrization of the observed light curve at low energies (discussed in the next section) 
and evaluated on t • P"' to take into account the delay due to a possible LIV effect, the 
measured signal spectrum Agjg and the effective area Agg. The PDF Pakg is composed of the 
uniform time distribution Afikg of the background events, the measured background spectrum 
Askg and the effective area Aes. No delay due to a possible LIV effect is expected in the 
background events of the ON data set. A(r„) and N' are the normalization factors of Pgjg 
and -Pskg respectively, in the [E, t) range of the likelihood £t. The coefficient Wg corresponds 
to the relative weight of the signal events in the total ON data set, derived from the number 
of events in the ON region tt-qn and the number of events in the OFF regions uoff weighted 
by the inverse number of OFF regions a. More details on the derivation of this function are 
given in Appendix IB. II 


(7) 

( 8 ) 
( 9 ) 
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3.3.2. Specific selection cuts and timing model 

The flare data set of the H.E.S.S. analysis (see section I2.ip was nsed with additional 
cnts. To perform the dispersion stndies, only nninterrnpted data have been kept. Thns, the 
analysis was condncted on the first 7 rnns, taken dnring the night of April 26th. Moreover, 
the cosmic ray flnx increases snbstantially for the 7th rnn, dne to a variation of the zenith 
angle dnring this night. This fact, along with its large statistical errors, leads ns to discard 
this rnn from the analysis. The 6th rnn shows little to no variability and was therefore 
also removed from the LIV analysis. Since within the ON data set, the signal and the 
backgronnd spectra have different indices (Tsig = 4.8 for the signal and TBkg = 2.5 for the 
backgronnd), the ratio S/B is expected to decrease with increasing energy. An npper energy 
cnt at Emax = 789 GeV was set, corresponding to the last bin with more than 3 a significance 
in the reconstrncted photon spectrnm (see the differential flnx dnring the flare in Fig. [1]). 
A lower cnt on the energy at Fmin = 300 GeV was nsed in order to avoid large systematic 
effects arising from high nncertainties on the H.E.S.S. effective area at lower energies. The 
intrinsic light cnrve of the flare, needed in the formnlation of the likelihood, can be obtained 
from a model of the timed emission or approximated from a snbset of the data. To be as 
model-independent as possible, it was here derived from a fit of the measnred light cnrve 
at low energies (with E < Ecut)- The high-energy events {E > Ecut) were processed in the 
calcnlation of the likelihood to search for potential dispersion. Here Ecut was set to Ecut = 
400 GeV, which is approximately the median energy of the ON event sample. Other cnts on 
the energy did not introdnce significant effects on the final resnlts. The histogram and the 
fit (Fig. [6]) were obtained as follows: the main idea was to preserve the maximnm detected 
variability in the PG 1553-1-113 flare, together with a significant response in each observed 
peak: 


• The binning was chosen so that at least two adjacent bins of the distribntion yield a 
minimnm of 3 a excess with respect to the average valne. 

• Simple parameterization have been tested on the whole data set (all energies): constant 
(x^/d.o.f=25/12), single Ganssian (x^/d.o.f=20/10) and donble Ganssian (x^/d.o.f=8.5/7) 
fnnctions. The latter is preferred, since it improves the qnality of the fit. This shape was 
chosen to fit the low energy snbset of events. Ghoosing a single Ganssian parametriza- 
tion wonld resnlt in a decrease of the sensitivity to time-lag measnrements by a factor 

of two. 

There is a gap of ~ 2 min between each two consecntive rnns. We did not consider the 
effect of these gaps as it is small with respect to the bin width of ~ 10 min. More importantly, 
their occnrrence is not correlated with the binning: one gap falls in the rising part of the 
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Time (s) 


Fig. 6.— Time distribution of the excess ON — aOFF in the hrst 6 runs (70971-70976), with 
energies between 300 GeV and 400 GeV. T = 0 corresponds to the time of the first detected 
event in run 70971. The vertical bars correspond to 1 cr statistical errors; the horizontal bars 
correspond to the bin width in time. The best fit, in red, was used as the template light 
curve in the maximum likelihood method; the ±1 a error envelope is shown in green. 
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light curve, one is at a maximum, two fall in the decreasing parts and none of the gaps is at 
the minimum. 

Table [H] in Appendix [R2] shows the number of ON and OFF events for the different cuts 
applied to the data. 


3.3.3. Results: limits on Tn and Eqq 

The maximum likelihood method was performed using high-energy events with Ei > E^ut- 
First, conhdence intervals (CIs) corresponding to 95 % conhdence level (1-sided) were de¬ 
termined from the likelihood curve at the values of where the curve reaches 2.71, which 
corresponds to the 90% C.L. quantile of a distribution. However, these CIs are derived 
from one realization only and do not take into account the “luckiness” factor of this mea¬ 
surement. To get statistically signihcant CIs (“calibrated CIs”), several sets were generated 
with Monte Carlo simulations, with the same statistical signihcance, light curve model and 
spectrum as the original data set. No intrinsic dispersion was artihciahy added. Each simu¬ 
lated data set produces a lower limit and an upper limit on r„. The calibrated lower (upper) 
limit of the conhdence interval is obtained from the mean of the distribution of the per-set 
individual lower (upper) limits. Both conhdence intervals (from the data only and from the 
simulated sets) are listed in Table [71 Sources of systematic errors include uncertainties on the 
light curve parameterization, the background contribution, the calculation of the ehective 
area, the energy resolution, and the determination of the photon index (see Appendix IB. 4|) . 

The resulting limits on the dispersion using the quadratic sum of the statistical errors 
from the simulations and the systematic errors determined from data and simulations were 
computed, leading to limits on the energy scale Eqq (Eq. [3]). The 95 % 1-sided lower limits 
for the subluminal case (s = -|-1) are: Eqg,i > 4.11 x 10^^ GeV and Eqq, 2 > 2.10 x 10^° GeV 

Table 5. Galibrated 95% 1-sided LL and UL (including systematic errors) on the 
dispersion parameter and derived 95% 1-sided lower limits on Eqq. 



Limits on 

Tn (sTeV ”) 

Lower limits 

on Eqg (GeV) 

n 

^^calib+syst 

^^calib+syst 

s= -1 

s= +1 

1 

-838.9 

576.4 

2.83 10^^ 

4.11 10^^ 

2 

-1570.5 

1012.4 

1.68 10^° 

2.10 10^° 
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for linear and qnadratic LIV effects, respectively. For the superInminal case (s = -1) the 
limits are: Eqg,i > 2.83 x 10^^ GeV and Eqg ,2 > 1-68 x 10^° GeV for linear and quadratic 
LIV effects, respectively. Fig. [7] shows a comparison of the different lower limits on Eqg,i 
and Eqg ,2 for the subluminal case (s = +1) obtained with AGN at different redshifts studied 
at very high energies. All these limits, including the present results, have been obtained 
under the assumption that no intrinsic delays between photons of different energies occur 
at the source. For the linear/subluminal case, the most constraining limit on Eqg with 
transient astrophysica l events has been obtained with GRB 090510: Eqg,i > 6.3 x 10^® GeV 
f Vasileiou et ah [2013). The mo s t con straining limits on Eqg with AGN so far have been 
obtained by Abramowski et ah ( 2011 ) with PKS 2155-304 data observed with H.E.S.S.: 
Eqg,i > 2.1 X 10^® GeV and Eqg, 2 > 6.4 x 10^® GeV for linear and quadratic LIV effects, 
respectively (95% GL, 1-sided). Gompared to the PKS 2155-304 limits, the limits on the 
linear dispersion for PG 1553-1-113 are one order of magnitude less constraining, but the 
limits on the quadratic dispersion are of the same order of magnitude since the source is 
located at a higher redshift. This highlights the interest in studying distant AGN, in spite 
of the difficulties due to limited photon statistics. 


4. Conclusions 

A VHE 7 -ray flaring event of PG 1553-1-113 has been detected with the H.E.S.S. 
telescopes, with a flux increasing of a factor of 3. No variability of the spectral index has 
been found in the data set, but indication of intra-night flux variability is reported in this 
work. In HE 7 rays, no counterpart of this event can be identihed, which may be interpreted 
as the sign of injection of high energy particles emitting predominantly in VHE 7 rays. Such 
particles might not be numerous enough to have a signihcant impact on the HE flux during 
either their acceleration or cooling phases. 

The data were used to constrain the redshift of the source using a new approach based on 
the absorption properties of the EBL imprinted in the spectrum of a distant source. Taking 
into account all the instrumental systematic uncertainties, the redshift of PG 1553-1-113 is 
determined as being z = 0.49 ± 0.04. 

Flares of variable sources can be used to probe LIV effects, manifesting themselves as an 
energy-dependent delay in the photon arrival time. A likelihood method, adapted to flares 
with a large amount of background and modest statistics, was presented. To demonstrate the 
analysis power of this method, it was applied to the H.E.S.S. data of a flare of PG 1553-1-113. 
This analysis relies on the indication of the intra-night variability of the flare at VHE. No 
signihcant dispersion was measured, and limits on the Eqg scale were derived, in a region of 
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Fig. 7.— Lower limits on Eqg,i from linear dispersion (left) and on Eqg ,2 from quadratic 
dispersion (right) for the subluminal case (s = +1) obtained with AGN as a function of 
redshift. Th e limits are given in terms of Epianck- The constraints fro m Mkn 421 have been 
obtain ed by Biller et ah 1 1999 1. from Mkn 501 by Albert et ah (2008), and from PKS 2155- 
304 by Abramowski et ah ( 2011 1. 
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redshift unexplored until now. Limits on the energy scale at which QG effects causing LIV 
may arise, derived in this work, are Eqg,i > 4.11 x 10^^ GeV and Eqg 2 > 2.10 x 10^° GeV for 
the subluminal case. Gompared with previous limits obtained with the PKS 2155-304 flare 
of 2006 July, the limits for PG 1553-1-113 for a linear dispersion are one order of magnitude 
less constraining while limits for a quadratic dispersion are of the same order of magnitude. 
With the new telescope placed at the center of the H.E.S.S. array that provides an energy 
threshold of several tens of GeV, a better picture of the variability patterns of AGN flares 
should be obtaine d. The future C herenkov Telescope Array (CTA) will increase the number 
of flare detections flSol et al.ll2013[l with better sensitivity, allowing for the extraction of even 
more constraining limits on the LIV effects. 
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A. Bayesian model used to constrain the redshift 

A Bayesian approach has been nsed to compnte the redshift valne of PG 1553+113 
in Section 13.21 The advantage of snch a model is that systematic nncertainties, which 
are important in Cherenkov astronomy, can easily be inclnded in the calcnlation. In the 
following, the notation 0 for the model parameters and + for the data set is adopted. 
All normalization constants are dropped in the development of the model, and the hnal 
probability is normalized at the end. 

Bayes’ Theorem, based on the conditional probability rnle, allows ns to write the poste¬ 
rior probability P(0|y) for the model parameters 0 as the prodnct of the likelihood P(+|0) 
and the prior probability P(0): 


P(0|F) oc P(0)P(F|0). 


The li kelihood is the qn antity that is maximized dnring determination of the best-£t 
spectrnm fjPiron et al.ll200ll) . It is at this step that the H.E.S.S. data, taken dnring the 
flare, were actnally nsed. The spectrnm model here is a simple power law corrected for the 
EBL absorption: 

(j) = Nx (E/Eo)-^ X 

The model parameters are then N, P and z. 


The prior is the most difficnlt and most interesting part of the model. To derive it, 
N and P are assnmed to be independent from each other and independent of the redshift. 
In contrast, the prior on the redshift might depend on N and P. Then, the prior can be 
simplihed nsing the conditional probability rnle: 


P(0) = P{z\N,T)P{N)P(T) 

As mnch as possible, weak assnmptions shonld be made to write a robnst prior then 
often flat (i.e. P oc const) are nsed. Priors shonld also be based on a physical meaning 
and not contradict the physical and observed properties of the objects. For the pnrpose of 
this model, the prior on N is assnmed to be flat and the prior on the spectral index is a 
trnncated Ganssian P(P) oc ^/^(r, Ppermi, ^r) if P < Ppermi and P(P) =oc const otherwise. 
The valnes of Ppermi and cxr are obtained by analyzing LAT data below 80 GeV (see section 
[3] and Table [3]). Here, it is assnmed that the intrinsic spectrnm in the VHE range cannot 
be harder than the Permz-LAT measnrement. ar takes into acconnt the statistical and 
systematic nncertainties on the Permz-LAT measnrement and also the systematic nncertainty 










29 


on the H.E.S.S. spectrum (a = 
(Jr = 0.33 for a mean value of F 


= 0.20, see lAharonian et al.ll2006bl ) added quadratically and 
Fermi 1.72. 


The prior on 2 ; is much more difficult to determine. A flat prior has no physical motiva¬ 
tions since the probability to detect sources at TeV energy decreases with the redshift. The 
number of sources detected at TeV energy is not sufficient to use the corresponding redshift 
distribution as a prior. 

A prior which takes into account the EBL, can be derived assuming a population of 
sources with a constant spatial density. In the small space element 47rz‘^dz, the number 
of such sources scales oc z‘^. For any given luminosity, their flux (which scales with the 
probability to detect them) is scaled by z~^exp{—T{z)). Lacking a proper knowledge of 
the intrinsic luminosity function of VHE 7 -ray blazars, a reasonable assumption on the 
detection probability of a blazar at any redshift is a scaling proportional to the flux for 
a given luminosity, i.e., oc z~‘^ exp{—T{z)). Putting everything together, the prior on the 
redshift reads P(z|V, T) = P{z) oc exp{—T{z)) 

Finally, the prior we use for our analysis is: 


F(0) oc exp(—r( 2 ;))A/G'(r, 1.72, 0.33) 


(Al) 


if T < 1.72 and 


P(0) oc exp{—T{z)) 


otherwise. Putting all the components of the model together and marginalizing over the 
nuisance parameters N and P, the probability on the redshift can be computed numerically. 
The obtained mean value is z = 0.49 ± 0.04. At a conhdence level of 95 %, the redshift is 
between 0.41 < z < 0.56. 


In this work, only the model of iFranceschini et al.l (120081) has been used. Other EBL 
models available in the literature predict slightly different absorption depths. This will lead 
to a small difference in the redshift. The use of a flat prior for the redshift distribution of the 
sources or a prior based on estimates of the HBLs luminosity function flAiello et al.l 120141) 
leads to changes of order of 0.01 on the resulting redshift. 
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B. Development of the LIV method 


B.l. Modified mciximum likelihood method 


In previous LIV studies with AGN flares flAlbert et al.l l2008l: lAbramowski et al.l 120111) 
the signal was clearly dominating over the background, whereas in the present study the 
signal-over-background ratio is about 2. The background has been included in the formula¬ 
tion of the probability density function (PDF): in the most general case, for given numbers 
of signal and background events s and b in the observation region (“ON” region), for a given 
dispersion parameter r„, the unbinned likelihood is: 


\ '^ON 

noFF\—yY\P{Ei,ti\s,b,Tn) (Bl) 
' 1=1 

The PDF P{Ei, ti\s, b, Tn) associated with each gamma-like particle characterized by its time 
ti and energy Ei contains two terms (signal and background): 


■ -Bsig(A'i , tj|Tfi)-|-(l Ws) ' Psksi^Eij ti) 


(B2) 


with 


Ws 


s 

s + b' 


(B3) 


rioN is the number of events detected in the source ON region included in the fit range 
[F^cut;-Emax] X [tmiriGmax]- 'H'OFF IS the number of events in the OFF regions, in the same 
{E, t) range; a is the inverse number of OFF regions. Pois(?7,oN|s 4- b) (Pois(?7,oFF|^/tt)) 
is the Poisson distribution with index tt-on (’^off) and parameter s + b [b/a). The likeli¬ 
hood function can be simplified by fixing s and b from a comparison of ON and OFF sets: 
s = n-oN — cmoFF and b = cmoFF- In this case, the Poisson terms in Eq. IB2l are equal to 1. 
The probabilities Psig and Pfilcg are defined as: 


Psisi.Piyti\Tji) • Psig(Pj, |t,^) 

\^n j 

PsksiPiPi) — ■ Psk^iEiyti) 


(B4) 

(B5) 


with 


Rs^g{E,t\Tn) 
RsksiP-i t) 


Etiue — 0 


^'true — 0 


D{E, Etme) AeffiEtrwi, t) Asig(£^tn.e) fsig(^ “ A. ' E'^^^^)dE 

E{E^ £/true) tj Aggg(i?(ggg) Aggg(i) 


(B6) 

(B7) 
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Pgig(£'j, tj|r„) is the probability that the event (Pj, ti) is a photon emitted at the sonrce and 
detected on Earth with a delay TnE'^. It takes into acconnt the emission (time distribntion 
Fgig(t) and energy spectrum Asig(P) at the source), the propagation (delay due to 

possible LIV effect) and the detection of a photon by the detector (H.E.S.S. energy resolution 
D{E, Etrne) and effective area Aef[{E, t)). -Pskg {Ei, ti) is the probability that the event {Ei, ti) 
is a background event; it is not expected to be variable with time, thus -^Bkg it) is a uniform 
time distribution: The background energy distribution Aakg is measured from 

OFF regions. N{Tn) (resp. N') is the normalization factor of the PDF Pgjg (resp. -Pfikg) in 
the range [Pcut; P ma v] x [tmin;tmax] where the likelihood fit is performed. 

Also, the energy resolution P(P, Ptrue) is assumed to be perfect in the range [Pcut; Pmax@- 
This leads to simplified expressions of Psig(Pj, fi|r„) and P^^^{EiAi)'- 

Psig(Pi,A|r„) = ] ■ Aeff(Pi,fi)Agig(Pi)Psig(fi - Tn ■ Pf) 

A [Tn) 

PsksiEiAi) ~ ■ Agff (Pj, fi)ABi^g(Pj)PBi^g 

The best estimate of the dispersion parameter r„ is obtained by maximizing the likelihood 

L(Tn). 


(B8) 

(B9) 


B.2. Selection cuts 

Table [6] shows the effect of the selection cuts on the number of ON and OFF events. 
Other choices of Pmin and Pcut did not introduce significant changes in the final results. 


®The actual energy resolution is of the order of 10 % in this range. 


Table 6. Selections applied to the ON and OFF data sets 


Selection 

# of ON events 

Weighted # of OFF events 

S/B 

Total sample 

461 (100 %) 

144.3 (100 %) 

2.2 

(1) = Time in 500-8500 s 

358 (77.7 %) 

95.8 (66.4 %) 

2.7 

(1) and P in 0.3-0.789 TeV 

154 (33.4 %) 

36.3 (25.1 %) 

3.2 

(1) and P in 0.3-0.4 TeV (Template) 

82 (17.8 %) 

14.2 (9.9 %) 

4.8 

(1) and P in 0.4-0.789 TeV (LH fit) 

72 (15.6 %) 

21.9 (15.2 %) 

2.3 








32 


B.3. Test of the method, confidence intervals 

The method has been tested on Monte Carlo (MC) simnlated sets. Each set was com¬ 
posed of rioN = 72 ON events, as in the real data sample: 

• s = 50 signal events with times following the template light curve (Fig. [6]) shifted by 
a factor r^^nj • E,; energies follow a power law spectrum of photon index Tsig = 4.8, 
degraded by the acceptance and convolved with the energy resolution. 

• b = 22 background events with times following a uniform distribution and energies 
drawn from a power law spectrum of index Enkg = 2.5, degraded by the acceptance 
and convoluted by the energy resolution. 

For a given injected dispersion, the maximum likelihood method is applied to each MC-simulated 
set. The initial light curve and energy spectrum were used as templates in the model instead 
of fitting them for each set. 

Figure [8] shows the means of the reconstructed dispersion versus the real (injected) 
dispersion for n = 1; for a given injected dispersion, error bars correspond to the RMS of the 
distribution of the best estimates ti. The blue line shows the result of a linear fit. The slope 
roughly corresponds to the percentage of signal in the total ON data set. It is due to the 
loss of sensitivity resulting from the part of the data sets with no dispersion. A systematic 
shift is observed of about 100 sTeV“^, well bellow la value - the RMS of the best estimate 
distribution is of 361 sTeV“^. The results in this paper have not been corrected for this 
bias. 

The coverage is not necessarily proper, i.e. the number of sets for which the injected 
dispersion value Tinj lies between the set’s lower limit (LL) and upper limit (UL) does not 
match the required 95 % 1-sided confidence level. The common cut used on the likelihood 
curves to get the LLs/ULs has been iteratively adjusted to ensure a correct statistical cov¬ 
erage: using this new cut, 95 % of the realizations provide CIs that include the injected 
dispersion r„^inj. The initial coverage was about 85 % for a cut on 21nT of 2.71. The new 
common cut, found iteratively at 3.5, ensures the desired 90 % 2-sided CL (approx. 95 % 
1-sided CL). Figure [9] shows the distributions of the best estimates, the 95% 1-sided LLs and 
ULs for Ti^inj = 0 sTeV“^ (linear case) and r 2 ,inj = 0 sTeV“^ (quadratic case); the means of 
the lower and upper limit distributions, shown as a blue vertical line, are used to construct 
the “calibrated confidence interval”. 

To get CIs from data, a maximum likelihood method is applied to the original data 
set and gives a best estimate The cut value determined from the simulations to 
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Fig. 8.— Means of the reconstructed dispersion versus the real (injected dispersion) for the 
linear case n = 1; for a given injected dispersion, errors bars correspond to the means of the 
distribution of the upper and lower limits (90 % 2-sided ~ 95 % 1-sided). The blue line is a 
linear fit to the points. The red line shows the ideally obtained curve Trecontmcted = 'Tinjected 
obtained in the case S/B = cx3. 

ensure proper coverage is applied on the original data set to obtain and The 

“calibrated” limits and combining from data together with MC results, 

are taken as 


r rcalib ^data i^MC j rMCi 

= T^best - l^hest - I 

/■/■rcalib _ data i | MC /■/■rMCl 

^ ^ — ^hest rbest ^ ^ I 


with and defined as the mean of the per-set best-estimate distribution, 

LL distribution, and UL distribution respectively. 


Table[7|lists the CIs determined in both ways, i.e., data-only and calibrated ones: 
and (resp. and are compatible within 10 %. In this work, calibrated 

CIs have been used to derive the final lower limits on Eqq. They are preferred over data-only 
CIs as they provides statistically well defined confidence levels. They al so ensure coherent 


compa rison with previous p ublished results , e.g. with PKS 2155-304 by lAbramowski et ah 
(2011) and CRB studies by Vasileiou et ah ( 2013 1. 
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B.4. Estimation of the systematics 

Estimations of the systematic effects on the dispersion measurement were performed. 
It was found that the main systematic errors are due to the uncertainties on the light 
curve parametrization. Other sources of systematic errors include the contribution of the 
background, effect of the change of photon index, the energy resolution and the effective area 
determination of the detector. To study the following four contributions, new simulated data 
sets have been built, each one with different input parameters: 

• background contribution: photons and background events have been reallocated within 
the ON data set in the fit range [Ecnt, .^max], introducing a la fluctuation in the number 
of signal event s in the ON data set; 

• effective area: set to a constant, equal to 120000 m^ for all energies and all times, 
which corresponds to a maximum shift of 10 % (the actual effective area increases with 
energy); 

• energy resolution: reconstructed energies have been replaced by the true energies; this 
corresponds to a shift of about 10 % on the reconstructed energy values; 

• photon index: changed by one standard deviation (±0.25). 

For the determination of systematic errors arising from the light curve parametrization, the 
calibration of the confidence intervals has been redone using successively the upper Icr and 
the lower la contours of the template, shown in Fig. El The change in mean lower and upper 
limits on the dispersion parameter gives an estimate of the systematic error associated 
to each contribution!^. An additional systematic contribution comes from the shift arising 
from the method found with simulation (see Appendix IB. 3 1) . Table [8] summarizes all studied 
systematic contributions. The overall estimated systematic error on is 330 sTeV“^ for 
the linear case (n = 1) and 555 sTeV“^ for the quadratic case (n = 2); they were included 
in the calculation of the limits on Eqg by adding the statistical and the systematic errors in 
quadrature. 


particular the errors on the peak positions constitute the most important part of the uncertainty 


on the template light curve contributing to the likelihood fit - see previous works, e.g. lAbramowski et al 


(|201ll) . Therefore, the covariance matrix of the fit of the template was studied in detail ; the peak positions 
were varied by values of ±lcr extracted from the covariance matrix. This study led to an increase in overall 
systematics of the order of 20% for ti and 40% for T 2 , and a decrease of maximum 7% and 2% of limits on 
Eqg,i and Eqg,2 respectively. 







Table 7. Linear (top) and quadratic (bottom) dispersion parameter; from left to right; 
best estimate, LL and UL from data (cut on likelihood curve), LL and UL from MC 
simulations (means of per-set LL and UL distributions), calibrated LL and UL 
(combination of data and MC), calibrated LL and UL including systematic errors. 
Dispersion parameters Tn^besu LLs and ULs are in sTeV"”. 


n 

^data 

n,best 



^MC 

n,best 


Ulmc 

^^calib 


with systematics 

1 

2 

-131.7 

-287.5 

-806.7 

-1449.9 

554.7 

853.6 

99.1 

217.2 

-526.3 

-942.0 

725.6 

1395.0 

-757.1 

-1446.7 

494.8 

890.3 

-838.9 

-1570.5 

576.4 

1012.4 


Table 8. Summary of all studied systematic contributions. The main systematic errors are 
due to the uncertainties on the light curve parametrization. 



Estimated error 
on input parameters 

Tl 

(sTeV-1) 

T2 

(sTeV-2) 

Background contribution 


< 45 

< 80 

Acceptance factors 

10 % 

< 1 

< 1 

Energy resolution 

10 % 

< 55 

< 85 

Photon index 

5% 

< 55 

< 50 

Light curve parametrization 


< 300 

< 500 

Systematic bias 


~ 100 

~ 200 

Total: 

< 330 

< 555 
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Fig. 9.— Distributions of the best estimates, the 95% 1-sided lower and upper limits from 
simulations in case of no injected dispersion (T^ inj = 0 sTeV”"'), for n = 1 (top) and n = 2 
(bottom); dispersion values are in sTeV“"^. The blue vertical line on the LL (resp. UL) 
distribution shows LL^'^ (resp. defined as the mean of the distribution. 
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